#if !defined(glib_g3_vector_h)
#define glib_g3_vector_h

#include "glib/operations.h"
#include <cmath>
//for IO operations #include<glib/g3/io.h>

namespace g3 {

template <class T>
class basic_vector;

typedef basic_vector<double> vector;
typedef basic_vector<float> vectorf;
typedef basic_vector<int> vectori;
typedef basic_vector<long double> vectorl;

template <class T>
class basic_vector {
public:
    T x, y, z;
public:
	static basic_vector<T> const nill; //(0,0,0) vector
public:
	typedef T * iterator;
	typedef const T * const_iterator;
public:
	T * begin()
		{ return &x; }
	const T * begin() const
		{ return &x; }
	T * end()
		{ return &z + 1; }
	const T * end() const
		{ return &z + 1; }
public:
	basic_vector()
		{ }
    basic_vector(T _x, T _y, T _z) 
		: x(_x), y(_y), z(_z) { }
	template <class U> 
	basic_vector(const basic_vector<U> & v) 
		: x(T(v.x)), y(T(v.y)), z(T(v.z)) { }
	basic_vector & zero() 
		{ x = y = z = 0; return * this; } 
public:
    basic_vector & operator += (const basic_vector<T> & v)
        { x += v.x; y += v.y; z += v.z; return * this; }
    basic_vector & operator -= (const basic_vector<T> & v)
        { x -= v.x; y -= v.y; z -= v.z; return * this; }
    basic_vector & operator *= (T a)
        { x *= a; y *= a; z *= a; return * this; }
    basic_vector & operator /= (T a)
        { x /= a; y /= a; z /= a; return * this; }
	// pagapov: the semantics is changed to lexicographical compare
	bool operator < (const basic_vector<T> & v) const;
	bool operator <= (const basic_vector<T> & v) const;
	bool operator == (const basic_vector<T> & v) const
		{ return x == v.x && y == v.y && z == v.z; }
	bool operator != (const basic_vector<T> & v) const
		{ return x != v.x || y != v.y || z != v.z; }
public:
	T & operator ()(int index)
		{ return * (&x + index); }
	const T & operator ()(int index) const
		{ return * (&x + index); }
	static int size() 
		{ return 3; }
public:
	/// returns a vector orthogonal to this
	basic_vector orthogonal() const;
	/// makes a non-zero vector of unit length
	basic_vector & normalize();
};

template <typename T>
bool basic_vector<T>::operator < (const basic_vector<T> & v) const
{
	if (x != v.x)
		return x < v.x;
	if (y != v.y)
		return y < v.y;
	if (z != v.z)
		return z < v.z;
	return false;
}

template <typename T>
bool basic_vector<T>::operator <= (const basic_vector<T> & v) const
{
	if (x != v.x)
		return x < v.x;
	if (y != v.y)
		return y < v.y;
	if (z != v.z)
		return z < v.z;
	return true;
}

template <typename T>
basic_vector<T> const basic_vector<T>::nill(0, 0, 0);

template <class T>
inline basic_vector<T> basic_vector<T>::orthogonal() const {
/*
	if (fabs(x) < fabs(y)) {
		if (fabs(x) < fabs(z)) {
			//x is the smallest component
			return basic_vector(y+z, -z-x, y-x);
		}
		else {
			//z is the smallest component
			return basic_vector(-y-z, x-z, y+x);
		}
	}
	else {
		if (fabs(y) < fabs(z)) {
			//y is the smallest component
			return basic_vector(-z-y, x+z, x-y);
		}
		else {
			//z is the smallest component
			return basic_vector(-y-z, x-z, y+x);
		}
	}
*/
	if (fabs(x) < fabs(y)) {
		if (fabs(x) < fabs(z)) {
			//x is the smallest component
			return basic_vector(0, -z, y);
		}
		else {
			//z is the smallest component
			return basic_vector(-y, x, 0);
		}
	}
	else {
		if (fabs(y) < fabs(z)) {
			//y is the smallest component
			return basic_vector(-z, 0, x);
		}
		else {
			//z is the smallest component
			return basic_vector(-y, x, 0);
		}
	}
}

template <class T>
inline basic_vector<T> & basic_vector<T>::normalize()
{
	T length = abs(* this);
	if (length != 0)
		* this /= length;
	return * this;
}


template <class T>
inline basic_vector<T> operator + (const basic_vector<T> & a, const basic_vector<T> & b) {
	basic_vector<T> t(a);
    return t += b;
}

template <class T>
inline basic_vector<T> operator - (const basic_vector<T> & a, const basic_vector<T> & b) {
	basic_vector<T> t(a);
    return t -= b;
}

//unary minus
template <typename T>
inline basic_vector<T> operator - (const basic_vector<T> & a) {
	return basic_vector<T>(-a.x, -a.y, -a.z);
}

template <class T>
inline basic_vector<T> operator * (const basic_vector<T> & v, T a) {
	basic_vector<T> t(v);
    return t *= a;
}

template <class T>
inline basic_vector<T> operator * (T a, const basic_vector<T> & v) {
	basic_vector<T> t(v);
    return t *= a;
}

template <class T>
inline basic_vector<T> operator / (const basic_vector<T> & v, T a) {
	basic_vector<T> t(v);
    return t /= a;
}

template <class T>
inline T scalar_product(const basic_vector<T> & v1, const basic_vector<T> & v2) 
	{ return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z; }

template <class T>
inline basic_vector<T> outer_product(const basic_vector<T> & v1, const basic_vector<T> & v2) {
	return basic_vector<T>(
		v1.y * v2.z - v1.z * v2.y,
		v1.z * v2.x - v1.x * v2.z,
		v1.x * v2.y - v1.y * v2.x
	);
}

//scalar triple product (v1,v2,v3) = (v1,[v2,v3])
template <class T>
inline T mixed_product(const basic_vector<T> & v1, const basic_vector<T> & v2, const basic_vector<T> & v3) 
	{ return scalar_product(v1, outer_product(v2, v3)); }

} // g3

template <class T>
class squarer<g3::basic_vector<T> > {
	T s;
public:
	explicit squarer(const g3::basic_vector<T> & v) : s(scalar_product(v, v)) { }
	operator T () const { return s; }
};

template <class T>
class absoluter<g3::basic_vector<T> > {
	T a;
public:
	explicit absoluter(const g3::basic_vector<T> & v) : a(::sqrt(sqr(v))) { }
	operator T () const { return a; }
};

#endif // glib_g3_vector_h
